The global distribution of South Africa’s four Anax species

There are four species of Anax (Emperor dragonflies) found in South Africa. Their ranges are extensive. This is a handy bit of code to visualize their distributions using R and GIS wizardry.

First, let’s set up our work space. These are the packages that we’ll need to run this code, check that they’re installed and in your library:

my_packages <- c("rmarkdown", "tidyverse", "sf", "rinat", "raster", "lwgeom", "leaflet","leafpop","rosm", "ggspatial", "prettymapr", "mapview", "rnaturalearth","readxl","htmltools","bslib") 
lapply(my_packages, require, character.only = TRUE)

Acquire Data

We’ll need some data. iNaturalist is a great place to find observational data on organisms around the world, including dragonflies. The data for this exercise was downloaded and stored in the GitHub repository - you can find it here.
The four species of dragonflies are:

Common Name Scientific Name
Blue Emperor Anax imperator Blue emperor
Orange Emperor Anax speratus Orange emperor
Vagrant Emperor Anax ephippiger Vagrant emperor
Black Emperor Anax tristis Black emperor


The species names will be used for the respective data sets below:
We’ll begin by reading in the data:

imperator <- read.csv("iNatData/obs_imperator.csv")
speratus <- read.csv("iNatData/obs_speratus.csv")
ephippiger <-read.csv("iNatData/obs_ephippiger.csv")
tristis <- read.csv("iNatData/obs_tristis.csv")


Then we’ll clean it up. Because iNaturalist makes use of citizen scientists to log their observations, the data isn’t always accurate or complete.
We’re filtering the data frames to remove observations without locations, and those with positional accuracy grater than 500m (dragonflies move around freely, so I’ve allowed for quite a bit inaccuracy here). We’ll also only include ‘research grade’ data, which has been somewhat verified by other users.

imperator <- imperator %>% filter(positional_accuracy < 500 & !is.na(latitude) & quality_grade == "research")
speratus <- speratus %>% filter(positional_accuracy < 500 & !is.na(latitude) & quality_grade == "research")
ephippiger <- ephippiger %>% filter(positional_accuracy < 500 & !is.na(latitude) & quality_grade == "research")
tristis <- tristis %>% filter(positional_accuracy < 500 & !is.na(latitude) & quality_grade == "research")


In order for the GIS wizardry to work, we’ll convert these tabular data frames to spatial objects using the Sf package:

imperatorSf <- st_as_sf(imperator, coords = c("longitude", "latitude"), crs = 4326) 
speratusSf <- st_as_sf(speratus, coords = c("longitude", "latitude"), crs = 4326) 
ephippigerSf <- st_as_sf(ephippiger, coords = c("longitude", "latitude"), crs = 4326) 
tristisSf <- st_as_sf(tristis, coords = c("longitude", "latitude"), crs = 4326)


Then we’ll use the following code to check that the data is indeed of the Sf class, and that the correct co-ordinate reference system is used for each:

cat("<details><summary>Click to expand the meaty results of these checks</summary>\n")
Click to expand the meaty results of these checks
class(imperatorSf) 

[1] “sf” “data.frame”

class(speratusSf) 

[1] “sf” “data.frame”

class(ephippigerSf) 

[1] “sf” “data.frame”

class(tristisSf)

[1] “sf” “data.frame”

st_crs(imperatorSf) 

Coordinate Reference System: User input: EPSG:4326 wkt: GEOGCRS[“WGS 84”, ENSEMBLE[“World Geodetic System 1984 ensemble”, MEMBER[“World Geodetic System 1984 (Transit)”], MEMBER[“World Geodetic System 1984 (G730)”], MEMBER[“World Geodetic System 1984 (G873)”], MEMBER[“World Geodetic System 1984 (G1150)”], MEMBER[“World Geodetic System 1984 (G1674)”], MEMBER[“World Geodetic System 1984 (G1762)”], MEMBER[“World Geodetic System 1984 (G2139)”], MEMBER[“World Geodetic System 1984 (G2296)”], ELLIPSOID[“WGS 84”,6378137,298.257223563, LENGTHUNIT[“metre”,1]], ENSEMBLEACCURACY[2.0]], PRIMEM[“Greenwich”,0, ANGLEUNIT[“degree”,0.0174532925199433]], CS[ellipsoidal,2], AXIS[“geodetic latitude (Lat)”,north, ORDER[1], ANGLEUNIT[“degree”,0.0174532925199433]], AXIS[“geodetic longitude (Lon)”,east, ORDER[2], ANGLEUNIT[“degree”,0.0174532925199433]], USAGE[ SCOPE[“Horizontal component of 3D system.”], AREA[“World.”], BBOX[-90,-180,90,180]], ID[“EPSG”,4326]]

st_crs(speratusSf) 

Coordinate Reference System: User input: EPSG:4326 wkt: GEOGCRS[“WGS 84”, ENSEMBLE[“World Geodetic System 1984 ensemble”, MEMBER[“World Geodetic System 1984 (Transit)”], MEMBER[“World Geodetic System 1984 (G730)”], MEMBER[“World Geodetic System 1984 (G873)”], MEMBER[“World Geodetic System 1984 (G1150)”], MEMBER[“World Geodetic System 1984 (G1674)”], MEMBER[“World Geodetic System 1984 (G1762)”], MEMBER[“World Geodetic System 1984 (G2139)”], MEMBER[“World Geodetic System 1984 (G2296)”], ELLIPSOID[“WGS 84”,6378137,298.257223563, LENGTHUNIT[“metre”,1]], ENSEMBLEACCURACY[2.0]], PRIMEM[“Greenwich”,0, ANGLEUNIT[“degree”,0.0174532925199433]], CS[ellipsoidal,2], AXIS[“geodetic latitude (Lat)”,north, ORDER[1], ANGLEUNIT[“degree”,0.0174532925199433]], AXIS[“geodetic longitude (Lon)”,east, ORDER[2], ANGLEUNIT[“degree”,0.0174532925199433]], USAGE[ SCOPE[“Horizontal component of 3D system.”], AREA[“World.”], BBOX[-90,-180,90,180]], ID[“EPSG”,4326]]

st_crs(ephippigerSf) 

Coordinate Reference System: User input: EPSG:4326 wkt: GEOGCRS[“WGS 84”, ENSEMBLE[“World Geodetic System 1984 ensemble”, MEMBER[“World Geodetic System 1984 (Transit)”], MEMBER[“World Geodetic System 1984 (G730)”], MEMBER[“World Geodetic System 1984 (G873)”], MEMBER[“World Geodetic System 1984 (G1150)”], MEMBER[“World Geodetic System 1984 (G1674)”], MEMBER[“World Geodetic System 1984 (G1762)”], MEMBER[“World Geodetic System 1984 (G2139)”], MEMBER[“World Geodetic System 1984 (G2296)”], ELLIPSOID[“WGS 84”,6378137,298.257223563, LENGTHUNIT[“metre”,1]], ENSEMBLEACCURACY[2.0]], PRIMEM[“Greenwich”,0, ANGLEUNIT[“degree”,0.0174532925199433]], CS[ellipsoidal,2], AXIS[“geodetic latitude (Lat)”,north, ORDER[1], ANGLEUNIT[“degree”,0.0174532925199433]], AXIS[“geodetic longitude (Lon)”,east, ORDER[2], ANGLEUNIT[“degree”,0.0174532925199433]], USAGE[ SCOPE[“Horizontal component of 3D system.”], AREA[“World.”], BBOX[-90,-180,90,180]], ID[“EPSG”,4326]]

st_crs(tristisSf)

Coordinate Reference System: User input: EPSG:4326 wkt: GEOGCRS[“WGS 84”, ENSEMBLE[“World Geodetic System 1984 ensemble”, MEMBER[“World Geodetic System 1984 (Transit)”], MEMBER[“World Geodetic System 1984 (G730)”], MEMBER[“World Geodetic System 1984 (G873)”], MEMBER[“World Geodetic System 1984 (G1150)”], MEMBER[“World Geodetic System 1984 (G1674)”], MEMBER[“World Geodetic System 1984 (G1762)”], MEMBER[“World Geodetic System 1984 (G2139)”], MEMBER[“World Geodetic System 1984 (G2296)”], ELLIPSOID[“WGS 84”,6378137,298.257223563, LENGTHUNIT[“metre”,1]], ENSEMBLEACCURACY[2.0]], PRIMEM[“Greenwich”,0, ANGLEUNIT[“degree”,0.0174532925199433]], CS[ellipsoidal,2], AXIS[“geodetic latitude (Lat)”,north, ORDER[1], ANGLEUNIT[“degree”,0.0174532925199433]], AXIS[“geodetic longitude (Lon)”,east, ORDER[2], ANGLEUNIT[“degree”,0.0174532925199433]], USAGE[ SCOPE[“Horizontal component of 3D system.”], AREA[“World.”], BBOX[-90,-180,90,180]], ID[“EPSG”,4326]]

cat("</details>\n")


Now that we’re confident that our data are accurate, research grade, and spatial objects, we can visualize them on a map:

ggplot() + 
  annotation_map_tile(type = "osm", progress = "none") + 
  geom_sf(data = imperatorSf, color = "#83aed1") + 
  ggtitle(expression("Distribution of " * italic("Anax") * italic(" imperator")))

ggplot() + 
  annotation_map_tile(type = "osm", progress = "none") + 
  geom_sf(data = speratusSf, color = "#e1a56d") + 
  ggtitle(expression("Distribution of " * italic("Anax") *  italic("speratus")))

ggplot() + 
  annotation_map_tile(type = "osm", progress = "none") + 
  geom_sf(data = ephippigerSf, color = "#afbd6c") + 
  ggtitle(expression("Distribution of " * italic("Anax") * italic(" ephippiger")))

ggplot() + 
  annotation_map_tile(type = "osm", progress = "none") + 
  geom_sf(data = tristisSf, color = "#9980ab") + 
  ggtitle(expression("Distribution of " * italic("Anax") * italic(" tristis")))


Great! These maps show the distribution of these species well, and the colour coding helps to differentiate them. However, these maps aren’t actually very useful to look at the observational data of one species in relation to another, nor does it allow us to visulise the data on any meaningful scale.
An interactive map, with all the species on it, would be much more useful! So we’ll use the Leaflet package to create a better map:

leaflet() %>%
  addTiles(group = "Default") %>%
  addCircleMarkers(data = imperatorSf, 
                   group = "Anax imperator", 
                   radius = 2, 
                   color = "#83aed1") %>%
  addCircleMarkers(data = speratusSf, 
                   group = "Anax speratus", 
                   radius = 2, 
                   color = "#e1a56d") %>%
  addCircleMarkers(data = ephippigerSf, 
                   group = "Anax ephippiger", 
                   radius = 2, 
                   color = "#afbd6c") %>%
  addCircleMarkers(data = tristisSf, 
                   group = "Anax tristis", 
                   radius = 2, 
                   color = "#9980ab") %>%
  addLegend(position = "topright", 
            colors = c("#83aed1", "#e1a56d", "#afbd6c", "#9980ab"), 
            labels = c("Blue Emperor", "Orange Emperor", "Vagrant Emperor", "Black Emperor"))


Much better!
Although, it would be easier to visualize this data if we could choose which species to display. So lets do that again, adding a layer control to choose which species’ data we’d like to see: And while we’re at it, it would make sense to see more details for each observation - perhaps even a clickable link to the iNaturalist page?
This chunk of code will create clickable links for us to add to popup labels:

limperatorSf <- imperatorSf %>% mutate(click_url = paste("<b><a href='", url, "'>Link to iNat observation</a></b>")) 
lsperatusSf <- speratusSf %>% mutate(click_url = paste("<b><a href='", url, "'>Link to iNat observation</a></b>")) 
lephippigerSf <- ephippigerSf %>% mutate(click_url = paste("<b><a href='", url, "'>Link to iNat observation</a></b>")) 
ltristisSf <- tristisSf %>% mutate(click_url = paste("<b><a href='", url, "'>Link to iNat observation</a></b>"))


And now we’re ready to create our user-friendly, and very pretty, interactive map:

leaflet() %>%
 addTiles(group = "Default") %>%
 addCircleMarkers(data = imperatorSf, 
                  group = "Blue Emperor", 
                  radius = 2, 
                  color = "#83aed1", 
                  popup = popupTable(limperatorSf,zcol = c("click_url"))) %>%
 addCircleMarkers(data = speratusSf, 
                  group = "Orange Emperor", 
                  radius = 2, 
                  color = "#e1a56d", 
                  popup = popupTable(lsperatusSf,zcol = c("click_url"))) %>%
 addCircleMarkers(data = ephippigerSf, 
                  group = "Vagrant Emperor", 
                  radius = 2, 
                  color = "#afbd6c", 
                  popup = popupTable(lephippigerSf,zcol = c("click_url"))) %>%
 addCircleMarkers(data = tristisSf, 
                  group = "Black Emperor", 
                  radius = 2, 
                  color = "#9980ab", 
                  popup = popupTable(ltristisSf,zcol = c("click_url"))) %>%
 addLayersControl( overlayGroups = c("Blue Emperor", "Orange Emperor", "Vagrant Emperor", "Black Emperor"), 
                   options = layersControlOptions(collapsed = FALSE) ) %>% 
 addLegend(position = "topright", 
           colors = c("#83aed1", "#e1a56d", "#afbd6c", "#9980ab"), 
           labels = c("Blue Emperor", "Orange Emperor", "Vagrant Emperor", "Black Emperor"))


It’s amazing that such small creatures have such a large range! There are at least 9 species of Anax dragonflies that are scientifically documented to migrate. While we still have much to learn about these impressive journeys, this distribution map shows that it’s likely that at least some populations of these four species migrate too.
A special mention to everyone who uses iNaturalist, citizen science makes observational data projects like this possible! And thanks to Jasper Slingsby for showing me how to make this, and for all his patience.